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Abstract. We present a detailed derivation of Fourier's law in a class of stochastic 
energy exchange systems that naturally characterize two-dimensional mechanical 
systems of locally confined particles in interaction. The stochastic systems consist of an 
array of energy variables which can be partially exchanged among nearest neighbours 
at variable rates. We provide two independent derivations of the thermal conductivity 
and prove this quantity is identical to the frequency of energy exchanges. The first 
derivation relies on the diffusion of the Helfand moment, which is determined solely by 
static averages. The second approach relies on a gradient expansion of the probability 
measure around a non-equilibrium stationary state. The linear part of the heat current 
is determined by local thermal equilibrium distributions which solve a Boltzmann-like 
equation. A numerical scheme is presented with computations of the conductivity 
along our two methods. The results are in excellent agreement with our theory. 
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1. Introduction 

A key to a comprehensive derivation of transport properties starting from a microscopic 
theory is to identify the conditions under which scales in a macroscopic system. This 
problem is especially challenging for interacting particle systems. To establish the 
conditions under which scales separate is to provide an understanding as to how a 
large number of particles whose microscopic motion is described by Newton's equations 
organise themselves in irreversible flow patterns at the macroscopic level. 

In the context of thermodynamics, macroscopic equations such as Fourier's heat 
law derive from conservation laws supplemented by phenomenological ones. The 
former concern say the conservation of mass or energy, and reflect the existence of the 
corresponding conservation laws at the level of Newton's equations. Phenomenological 
laws on the other hand provide linear relations between the currents associated with 
the flow of conserved quantities and thermodynamic forces in the form of gradients of 
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the conserved quantities, thereby introducing a set of transport coefficients. Though 
the values of these coefficients can be precisely measured and are usually tabulated 
for the sake of their use in the framework of applied thermodynamics, they cannot be 
otherwise determined without knowledge of the underlying dynamics. A ffist-principles 
based computation of the transport coefficients consequently requires a deep knowledge 
and understanding of the dynamics, as well as their statistical properties. 

For this purpose, a common procedure in non-equilibrium statistical physics -see 
e.g. Uhlenbeck's discussion of Bogoliubov's approach to the description of a gas in [3] 
or van Kampen's views on the role of stochastic processes in physics [1]- is to apply the 
two-step programme which consists of (i) identifying an intermediate level of description 
-a mesoscopic scale- where the Newtonian dynamics can be consistently approximated 
by a set of stochastic equations, and (ii) subsequently analyzing the statistical properties 
of this stochastic system so as to compute its transport properties. 

The ffist part of this programme was successfully completed in a recent set of papers 
[H E], in which we introduced a class of Hamiltonian dynamical systems describing 
the two-dimensional motion of locally confined hard-disc particles undergoing elastic 
collisions with each other. Whereas the confining mechanism prevents any mass 
transport in these systems, energy exchanges can take place through binary collisions 
among particles belonging to neighbouring cells. Under the assumption that binary 
collisions are rare compared to wall collisions, it was argued that, as a consequence 
of the rapid decay of statistical correlations of the confining dynamics, the global 
multi-particle probability distribution of the system typically reaches local equilibrium 
distributions at the kinetic energy of each individual particle before energy exchanges 
proceed. This mechanism naturally yields a stochastic description of the time evolution 
of the probability distribution of the local energies in terms of a master equation to 
be described below. An important property in that respect is that the accuracy of 
this stochastic reduction can be controlled to arbitrary precision by simply tuning the 
system's parameters. 

In [U [2] , arguments were presented supporting the result that the heat conductivity 
associated with this master equation has a simple analytical form, given by the frequency 
of energy exchanges. In a sense, this result is similar to the fact that, for instance, 
uniform random walks describing tracer dynamics have diffusion coefficients given in 
terms of jump probabilities of the walkers. However, contrary to such systems, the 
stochastic system at hand lacks a special property which facilitates the derivation of 
such results, namely the gradient condition. A system obeys the gradient condition 
if there exists a local function such that the current, whether of mass or energy, can 
be written as the difference of this function evaluated at separate coordinates [5]. In 
such cases, the diffusion coefficient is determined through a static average only, and is 
therefore easy to compute. Given that our stochastic system does not verify the gradient 
condition, it is in fact remarkable that the heat conductivity should take such a simple 
form. 

In this paper, we complete the second part of the programme described above and 
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provide a systematic derivation of the heat conductivity associated with heat transport 
in the class of stochastic systems derived in [H [2]. We justify in particular why, in spite 
of the fact that our systems do not obey the gradient condition, the heat conductivity 
can still be computed from the Green-Kubo formula through a static average only. We 
achieve this by an alternative method which consists in considering the stationary heat 
flux produced by a temperature gradient across the system. We suppose the system 
is a two-dimensional slab. Along the first dimension, the system has finite extension, 
with the corresponding borders in contact with stochastic thermal baths at different 
temperatures. The system may be taken to periodic along the second dimension. As a 
result, a temperature gradient develops along the first dimension, which induces a heat 
flux from the warmer border to the colder one. We then set up a scheme to compute 
the resulting non-equilibrium stationary state and obtain the linear relation connecting 
the temperature gradient to the heat flux. This scheme proceeds by consistently 
solving the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy up to pair 
distributions through a type of Chapman-Enskog gradient expansion. As the heat 
current associated with our systems depends on neighbouring pairs only, the knowledge 
of the pair distribution function to first order in the temperature gradient makes possible 
the computation of the heat current in the linear regime and, consequently, that of the 
thermal conductivity, thereby completing the derivation of Fourier's law. 

Our theoretical results are furthermore supported by the results of direct numerical 
simulations of the stochastic system, with excellent agreement. Simulation methods of 
kinetic processes governed by a master equation have been extensively developed after 
Gillespie's original work [HI [7]. Here we describe a method which accounts for continuous 
energy exchanges among all the pairs of neighbouring cells. 

The paper is organised as follows. The master equation is described in section 
[21 with relevant definitions. Statistical objects are described in section [3], identifying 
equilibrium and non-equilibrium stationary states. Local temperatures are defined in 
section H] and the energy conservation law derived, thereby identifying heat currents. 
Section [5] offers a first computation of the heat conductivity through the method of 
Helfand moments. An alternative computation is presented in section [6l solving the 
BBGKY hierarchy as sketched above. The numerical scheme for the computation of the 
heat conductivity is discussed in section [71 with a presentation of the results. Finally, 
conclusions are drawn in section [H] 

2. Master equation 

Consider a lattice of confining two-dimensional cells, each containing a single hard-disc 
particle, and such that particles in neighbouring cells may perform elastic collisions 
with each other, thereby exchanging energy. Such mechanical systems with hard-core 
confining mechanisms have recently been considered in [U [2] (see figure [I]). In these 
systems, one distinguishes the local dynamics from the interacting dynamics. On the 
one hand, the local dynamics are characterized by a wall-collision frequency, z/w, which 
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depends on the geometry of the confining cell as well as on the kinetic energy of the 
moving particles. On the other hand, the interacting dynamics are characterized by 
the frequency of binary collisions, i.e. the frequency of collisions between neighbouring 
particles, z/b- 




Figure 1. Snapshot showing an example of an interacting particle system with a 
confining mechanism, such as described in [2 . The coloured particles, colour-coded 
according to their kinetic energies from blue to red, move among an array of fixed black 
discs. The system has two parameters, the diameters of the fixed and moving discs. 
They are chosen so that their sum is larger than the distance between neighbouring 
fixed discs. This guarantees that the moving discs are trapped within their cells; their 
centers move within semi-dispersing billiards of bounded horizon (exterior intersection 
of the black circles). Though they are confined, the moving discs can nevertheless 
exchange energy through binary collisions as long as their diameters are properly 
chosen. The numbers indicate the values of the kinetic energies in every cell. In 
the limit where binary collisions are rare with respect to wall collisions, the energy 
exchange dynamics reduces to the stochastic evolution described by equation ([1]). 

Under specific conditions, the scale separation, z/b -C z/w, is achieved, which is 
to say that individual particles typically perform many collisions with the walls of 
their confining cells, rattling about their cages at higher frequency than that of binary 
collisions. In this regime, Liouville's equation governing the time evolution of phase- 
space densities reduces to a master equation for the time evolution of local energies. 
Moreover the validity of this reduction is controlled by the scale separation between 
the two collision frequencies z^b and z/w, and becomes exact in the limit of vanishing 
binary collision frequency. On the contrary, in the absence of a clear separation between 
binary and wall collision frequencies, the dynamics is not reducible to such a stochastic 
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description and cannot be addressed in the context of this paper -see [2] for a discussion 
of that hmit. 

Throughout this article, we assume the vahdity of the scale separation between the 
wall and binary collision frequencies and focus on the stochastic reduction of the energy 
exchange dynamics, considering, as our starting point, the mesoscopic level description 
of the time evolution of probability densities as a stochastic evolution. The time 
evolution is thus specified by a master equation which accounts for the energy exchanges 
between neighbouring cells. Note that the dimensionality of the dynamics specifies the 
maximum dimension of the energy cell array we may consider. In the case of the 
dynamics depicted in figure [U, the array of energy variables would be two-dimensional. 
However, without loss of generality, we can simplify the description and consider instead 
a one- dimensional array according to which every cell has two neighbours instead of four, 
one on each sid^. A similar construction can be carried out in three dimensions, starting 
with systems of confined hard balls. 

We can thus leave aside the underlying dynamics and consider instead a system of 
cells along a one-dimensional axis with energies {ei, . . . , e^} and let Pn{^i, ■ ■ ■ , ^n, t) 
denote the time-dependent energy distribution associated with this system. The time 
evolution of this object is specified according to 

dtPN{eu...,eN,t) = (1) 

)Pn{. . . ,ea + rj,ea+i -rj,...,t) 

a=l •' 

-W{ea,ea+i\ea - r],ea+i + r])PN{. . . ,ea,ea+i, ■ ■ ■ ,t) , 

where, for the sake of the argument, we assume periodic boundary conditions and 
identify cells a = + 1 and a = 1. 

The kernel W specifies the energy transition rates, whereby a pair of energies {ea, e^} 
exchange an amount rj of energy. In the case of systems such as shown in figure [H it is 
given by 

W{ea, eb\ea -r],eb + v)= \2\r 7^ / ^^P^^ / dvadvb (2) 

X eab ■ Vab 5 (^ea - y S (^Eb - ^^6^ ^ (v - -^[{^ab " - (Bab " Vb)^]j , 

where Cab denotes the unit vector joining particles a and b with respective velocities Va 
and Vb, 4> is the angle between the direction of this unit vector and a reference axis, 
R denotes the position of the center of mass of the two particles, m their masses, pm 
their radii, and |£p_p„^(2)| is the configuration-space volume which they occupy, with p 
a parameter characterizing the geometry of the confining cell. Thus t] is the amount of 

I Though it is of course possible to consider a two-dimensional array, as for instance shown in figure 
[1] we will focus our attention on onc-dimensional arrays because, in the presence of a temperature 
gradient along, say, the horizontal direction, no conduction takes place in the vertical direction. The 
only relevant energy exchange processes happen along the direction of the temperature gradient. 
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energy exchanged by the two particles of respective energies and ef, in the colhsion 
process. 

The expression ([2]) of extends beyond the mechanical systems described above. 
It applies to all mechanical models in which locally confined particles interact with their 
nearest neighbours through hard-core collisions. The dimensionality of the dynamics, 
whether two or three, is of course relevant to the specific form of the kernel, in particular 
as far the computation of the velocity integrals goes. However, whether the underlying 
dynamics is two- or three-dimensional, it can be shown that the transport properties of 
the energy exchange process to be established below are obtained in similar ways. 

In the above expression, one shows that the spatial integral / d(f)dR decouples 
from the velocity integral /- ^.^ ^^^dVadvi,. The former quantity represents the volume 
that the center of mass of particles a and b can occupy given that the two particles 
are in contact, performing a collision. Now rescaling the time variable by a factor 
2pjn/ d(f)dR/[^y7lm\Cp^p^{2)\], which amounts to converting time to the units of the 
inverse of the square root of an energy, we can rid equation ([2]) of all geometric factors 
and write the expression of the kernel in terms of Jacobi elliptic functions of the first 
kind, denoted K. Assuming ea < e?,, the kernel takes the expression 

Wiea,e,\ea-v,eb + v) = d:^y< \ ^k{^) , - e, < r/ < , (3) 

iv^^(^)' 0<,<e.. 
For efo < ea, the kernel is defined in a similar way, using the symmetry of W with respect 
to its arguments, i.e. exchanging and and changing the sign of rj. 

The rate at which two neighbouring cells with energies ea and ef, exchange energy is 
obtained by integrating W over the range of possible values of the amount rj of energy 
exchanged, 

u{ea,eb)= dr]W{ea,eb\ea - 'ni(^b + v) , (4) 

which, using equation ([3]), is easily computed in terms of first and second kind elliptic 
functions, K and E, to read 




2E[- \ - [l~-\K [- 



eh J \ e^J \eb 



(eb > ea) • (5) 



This is a symmetric function of ea, e?,, i.e. i/(ea, Cb) = z/(efc,ea) if ea > eh- It is 
homogeneous in the sense that, given a positive constant a, we have z/(aea,aeb) = 
A/az/(ea,efe). 

Likewise, the average heat exchanged between two cells at respective energies 
and eb, is 

j{ea,eb)= / dr]r]W{ea,eb\ea-ri,eb + ri) , 




(6) 
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which is an antisymmetric function, i.e. defined by j{ea, Cf,) = —j{eb, e^) when > e?,. 
3. Stationary states 

A stationary state of ([1]) is a time- independent probabihty distribution Pj^^\ei, . . . , e^) 
such that, for all pairs {ea,eb), 

J dr] W{ea + r],eb- r]\ea, eb)PN''\- + ?7, e;, - ?7, .. .) 

-Wiea, eb\ea - r^, e, + r])Pi^'^ (..., e„, e,, ...)]= . (7) 

3.1. Equilibrium states 

Equilibrium measures have densities which depend on the local energies through their 
sum only, e.g. the micro-canonical distribution 6{ei + . . . + ej\f — E) and the canonical 
distribution exp[—j3{ei + . . . + e^)]. The former is associated with an isolated A^- 
cells system, the latter to a system in contact with a thermal bath at fixed inverse 
temperature (3 = 1/T. Note that we assume the temperature is measured in the units 
of the energy, or, equivalently, that the Boltzmann constant is unity, k-Q = 1. 

The equilibrium average of the rate of energy exchanges (jlj) is the collision frequency, 
which in the canonical ensemble at temperature T, is given by 
1 f 

vb = j deade5Z/(ea,eb)exp[-(ea + eb)/T] = Vt. (8) 

Hence the choice of the time scale. 

We further notice that, for finite size systems of cells and total energy E = NT, 
i.e. average energy T per cell, the collision frequency is obtained from the micro- 
canonical average, 

1 pNT rNT-ea / f A_(:\N-3 

"-=(mL ""L de..fe.aA'-l)(A'-2)(l-^) . 

= Vt + o{N'^). (9) 

For future sake, we also notice the identity 

J dxdy{x — y)'^iy{x, y) exp(— x — y) = 3 . (10) 

3.2. Non- equilibrium states 

A non-equilibrium stationary state occurs e.g. when the two ends of a one-dimensional 
channel are put in contact with thermal baths at inverse temperatures /3_ and with 
7^ so that a stationary heat fiux will fiow from the hot to the cold end. One might 
hope that the corresponding stationary state would have the local product structure 

Pj'^(ei, . . . ,eAr) = /3i ■ ■■ f3Nexp[-(3iei - ... - /?AreAr] , (11) 

with a temperature profile specified according to Fourier's law dx[n{x)dxT{x)] = 0. 
However such a distribution is not stationary since an energy exchange displaces an 
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amount of energy between two cells at different local temperatures. We will however 
see in Section [6] that (fTTl) is indeed a good approximation to the actual non-equilibrium 
stationary state and can be used to compute the linear response. 



4. Conservation of energy 

In order to analyze heat transport in our system, we need to consider the time evolution 
of local temperatures 

kBT{a,t) = (ea) . (12) 

We have, using equation ([1]), and assuming we have a one- dimensional lattice of cells so 
that only terms (a — 1, a) and (a, a + 1) contribute, 

dtkBT{a, t) = j dei . . . deAre„9tPiv(ei, . . . , e^v, t) , 

= \Y^ /dei . . . deNd7]\eaW{ea' + r/, e?,/ - 7]\ea' ,eb')PN{- . . ,e'^ + r], . . . ,eb' - r], . . . ,t) 



a',b' 



I 



dei 



t) 



dei 



(13) 



-eaW{ea',eb'\ea' - r],eb' + r])PN{- ■ • ,ea', • • • . . . ,t) , 
deNdrj eaW{ea + r], )Pn{- • • ,ea + ?7,ea+i - V, 

-eaW{ea, ea+i\ea - V, ea+1 + r])PNi- ■■,ea, ea+i, ...,t) 

-(^aW{ea-i, ea\ea-i -r],ea + v)Pn{- ■ ■ , Ca-i, ea,...,t) , 
deArd?7 [ - T]W{ea, ea+i |e„ - i], e^+i + r])PNi- Ca+i, ■■■,t) 

+riW{ea-i, ea\ea-i -r],ea + 'n)PNi- ■ ■ , Ca-i, ea,...,t) . 

The two terms on the RHS of the last line denote the energy currents flowing between 
the neighbouring cells, 

dtkBT{a,t) = Ja-l,a{t) - Ja,a+l{t) , (14) 

where 

Ja,b{^) = J dei... de^jiea, eb)PAr(ei, . . . , ejv, t) (15) 

denotes the average energy flux from cell a to cell b with respect to the distribution Pm- 
Thus equation (1141) is a conservation equation for the energy. We remark here that we 
have dropped the dimensional factors, setting the cells lengths to unity. 

Fourier's law relates the heat current to the local gradient of temperature through 
a linear law which involves the coefficient of heat conductivity. We now set out to derive 
Fourier's law and thereby compute this quantity. 



5. Helfand moment 



The heat conductivity which is associated with heat transport can be computed by 
considering the linear growth in time of the mean squared Helfand moment for that 
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process, which undergoes a deterministic diffusion in phase space |H]. This is a 
generahsation of Einstein's formula which relates the diffusion coefficient to the mean 
squared displacement. 

Assuming a one-dimensional lattice of cells with unit distance between neighbouring 
cells, we may write the Helfand moment associated with heat transport as 

N 

H{t) = Y^aea{t). (16) 

a=l 

Changes in this quantity occur whenever an energy exchange takes place between any 
two neighbouring cells. Thus considering a stochastic realisation of the energy exchange 
process, we have a sequence of times r„ at which a given pair of cells A;„ + 1) 
performs an energy exchange ??(efc„, efe„+i), i.e. ek„ ek„ - i] and ek„+i efc„+i + r/. 
The corresponding Helfand moment therefore evolves according to 

H{Tn) = H{Tn-i) + r/(efe„, efc„+i) , (17) 

and has overall displacement 

n 

AH{t,,) = HiT„) - H{t,) = 5]r/(6fc,,e,,+i) . (18) 

i=l 

The thermal conductivity can be computed in terms of the equilibrium average of 
the mean squared displacement of the Helfand moment according to 

K = lim , \ lim (—AH(Tn?) , (19) 

where {.)e/n denotes a micro-canonical equilibrium average of the cells system at 
energy E = NT, which involves both energy exchanges and relaxation times between 
then. When is large, boils down to an average with respect to the canonical 
distribution at temperature T. In equation (fT9l) . the time r„ corresponding to n events 
is, by the law of large numbers, the typical time it takes for n collision events to occur 
in a system of A^ cells, which is easily expressed in terms of the collision frequency (IHl), 

Ti 

lim — = Nub . (20) 

' n 

Substituting AH from equation (fTSl) into equation (fT9l) . we obtain the thermal 
conductivity. 



K = lim lim 

NT^ 



5:(-^r/(6,^,e,,+i)n (21) 



2r 

n—l n I Y 



+2E E (^^(efc,,efc,+i)r7(efc^,,efc^,+i; 



i=l jr=i+l 

which is equivalent to the corresponding result using the Green-Kubo formula. 

The first of the two terms in the brackets on the RHS of ( !2TI) is a sum of n identical 
static averages, which, for large A^, is approximated by the canonical equilibrium average 

r7(ea,efe)^^^ = — ^ ^ deade^ d?7?7^W^(ea, e^lca - r^) exp[-(ea eb)/r] , 

= 2T2 . (22) 
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Notice here that, in order to compute the second moment of the average energy 
exchanged between the two cells, we divided W hj u^, the rate at which these exchanges 
occur. 

The second term on the RHS of equation (12T!) . on the other hand, is a sum of 
dynamic averages which are generally difficult to compute. We contend that it goes 
to zero when ^ oo. This is to say that cross- correlations of energy transfers 
{ri{ek^, efc^+i)?7(efc^., efc^.+i))^ exist only so long as the system size is finite. The reason for 
this is that, in large enough systems, new energies keep entering the dynamic averages 
in equation (I2T!) . as if the systems were in contact with stochastic reservoirs, which is 
enough to destroy these correlations. 

In the following section, we turn to the evaluation of the non-equilibrium stationary 
state and will provide more definitive arguments that concur with these heuristic 
reasoning. In section [7] we discuss the results of numerical computations of the Helfand 
moment which provide further insight into the decay of these dynamic averages as 
iV ^ oo. 

Accepting the claim that dynamic averages vanish in equation (1211) and thus 
that only static averages contribute to the heat conductivity, we conclude that the 
thermal conductivity associated with the process defined by ([1]) is equal to the collision 
frequenc}|§|, which, for a general equilibrium temperature T, is 

k = ub = Vt. (23) 
6. Chapman-Enskog gradient expansion 

Consider the marginals of the A^-cell distribution function, given by the one- and two-cell 
distribution functions, 

pW(e„,t) = J nde,Pjv(ei,...,...,ejv,t), 

The time evolution of the former can be written in terms of the latter, 
9tPa^\e,t) = J Y[deidtPN{ei,---,ea-i,e,ea+i,...,eN,t), 

= J de'dv[wie + r],e'- r]\e, e')Pi,li(e + r^, e' - r^) 

-Wie,e'\e-v,e' + v)Pi%ie,e') 

+ W{e + v,e'- v\e, e)P^L{e + r^, e' - r^) 

-W{e,e'\e-v,e' + v)Pl'l,{e,e')]. (25) 

This equation is the first step of the BBGKY hierarchy. We want to solve this equation 
in the non-equilibrium stationary state which is generated by putting the system 

§ The dimensions of these quantities differ by a length squared -the separation between the energy 
cells- which we have set to unity. 



(24) 
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boundaries in contact with thermal reservoirs at two different temperatures, T_ 7^ T+, 
with respective distributions Pl{^) = /3_ exp(— /3_e) and -PR(e) = /3+ exp(— /3+e), at the 
corresponding inverse temperatures /5_ and 

In the presence of such a temperature gradient, we expect that a uniform stationary 
heat current, as defined by equation (fTSll and henceforth denoted Jh, will establish itself 
throughout the system. Provided the local temperature gradient (T+ — TJ)/N is small, 
this current would be given, according to Fourier's law, by the product of the local 
temperature gradient and heat conductivity n\ 

Jn = {Ta+i - Ta) + 0[(T+ - T_)/iV]2 . (26) 

This equality should hold for any 0<a<N,a = and a = N + 1 corresponding to 
the two thermal reservoirs. Notice that the heat conductivity is itself a function of the 
temperature and thus varies from site to site. 

In order to derive Fourier's law ( l26i) . as well as the expression of the heat 
conductivity k, and thus establish equation (l23l) . we compute the heat current using 
equation (fTSil . which, in terms of the two-cell distribution P„ Q^_|_^(e, e'), becomes 



Jn^J dede'j(e, e')Pi,li(e, e') . (27) 
As noticed above, Jh on the LHS of this equation is defined to first order in the local 

(2) 

temperature gradients. We therefore need to compute PaJ+i{^, e') to that order as well, 
which is to say it must be a stationary state of equation ( l25i) up to second order in the 
local temperature gradient. 

In order to find this approximate stationary state, we perform a cluster expansion 
of the two cell distribution function and start by assuming that it factorises into the 
product of two one-cell distribution functions, 

i^i?(e,6')=Pi^H6)P«(6'). (28) 

This solution is as yet incomplete and must be understood as being only the first step 
in obtaining the consistent solution which is to be written down below. 

Substituting equation ( |28l) into equation ( |25l) . it thus reduces to a Boltzmann-Kac 
equation, which, owing to the symmetries of the kernel, can be written as 

9,P«(e,t) = I d6'dr^W^(e,e'|e-r/,e' + r^)[P«(6 + r^)P«i(e'-r^) (29) 

A first approximation to the stationary state of this equation is given by the local 
thermal equilibrium solution, 

P«(e) = /5„e-^-. (30) 

By definition of the local temperature, we have 

1 

deeP«(e)=/?-i = T,. (31) 
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One verifies that the form fl30l) is an approximate stationary solution of the Boltzmann 
Kac equation fl29|) . Indeed, let Pa = P denote the equilibrium temperature, Pa±i = 
f3 ± 6(3. We have 



f3''e-^^'-'''>7]6(3 + 0{6(3''), (32) 
Likewise 

p!.'\^ + v)P^'Me -V)- Pl'\e)Pi'M^') = -(5\-^^^^'\6(5 + 0{6(?) .(33) 

Therefore, to 0(5/3^), the sum of equations ( l32l) and ( l33l) vanishes. This is to say the 
local thermal equilibrium distributions (!30|) give approximate stationary states of the 
Boltzmann-Kac equation (l29l) to that order. 

Going back to the two-cell distribution, equation fl28l) . we proceed to compute the 
next order of the BBGKY hierarchy. One might hope that local equilibrium solutions 
are also a good approximation to the two-cell distribution function. This is however not 
the case. We can see this by considering the time evolution of the two-cell distribution. 



dtP, 



(2) 



a,a+l V-) 



e, e 



(34) 



dr/W^(e, e'|e - r^, e' + ri) [pfl^^{t + r^,e'-r^)- P^l,{e, e') 
+ j de"dvW{e, e"\e - r/, e" + [pS,,,,+i(e" -v,e + v, e') - PS,,,,+i(e", 6, e') 

+ J de"dvW{e\ e'y - r/, e" + v) [Pi2+i,a+2(e, e' + r^, e" - r/) - Pi,li,,+2(e, e', e" 

We momentarily suppose that the three-cell distributions can be consistently written as 
a product of one-cell distributions, 

p!SA^. e") = P,«(6)Pi^^(e')P,«(e") . (35) 

Plugging this form into equation ( l34l) and assuming a stationary state, we arrive, after 
expanding all terms to 0{6f3), to the equation 

= /3''e-^^'+''^6/3(^jie, e') - /3 J de"e-^^"j(e, e") + P J de" e-^'" ]{e' , e")} • (36) 

This is however a contradiction since the current does not have the gradient form which 
(!36|) implies. 



j(e, e')^(3j de"e~^''j{e, e") - f3 J de"e-^'' j{e' , e") . 



(37) 



Indeed, the only possible solution of this equation is a current given in terms of the 
difference of two local functions, say j(e, e') = /(e) — /(e'), which, clearly, equation ([6]) 
does not satisfy. 

To work around this difficulty, we must go back to equation (l28l) and perform a 
cluster expansion to include 0{6j3) corrections in the form 

e') = Pi^)(e)Pii\(6') + (/?,±i - /?a)QHi(e, e') , 



lTW±^9S±i(/5e,/?e') 



(38) 
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where the first two terms in the second line come from the expansion of the one-cell 
distributions and g*-^-* is one plus an expression derived from Q^'^\ 

Notice that, for definiteness, we require that / de'P^^^^-* (e, e') = Pj^^(e), implying 

/de'gS(e,e')=0. (39) 

Substituting the form ( !38|) into the RHS of equation ( !25|) . we must have the 
cancellation of all 0{Sl3) terms, and thus 

= J de'dr7e-^^V(e,e'|e-r/,e' + r7){gi5+i[/3(e + r/),/3(e'-r/)] (40) 
-q^a}+im (3e') - gS-i[/3(e + r^) ,/3(e' - v)] + gH-i(/5e, /^e')} . 
Given that the system is large, we expect a-ti to converge to the forms 

(la,a+i{x, y) = q{x, y) + 0{5I3) , ^^^^ 
£-1(2;, y) = q{y, x) + 0{6f3) . 

If so, let us consider q to be the sum of a symmetric function g and an antisymmetric 
function h, q{x, y) = g{x, y) + h{x, y), with g{x, y) = g{y, x) and h{x, y) = —h{y, x). It 
is obvious that the terms between the brackets of equation (HOl) cancel each other with 
respect to the symmetric part g. As far as the antisymmetric part, on the other hand, 
we must have 



dye' 



. (42) 



dzW{x, y\x — z,y + z)h{x + z,y — z) — z/(x, y)h{x, y) 

We argue that the only possible solution of equation (l42l) is a function h which 
depends on the sum of its arguments, which is obviously not an antisymmetric function. 
Therefore g is a symmetric function of its arguments, 

q{y,x) = q{x,y) . (43) 

With the two-cell distribution fl38|) . we write the three-cell distribution, now 
including the second order terms of the cluster expansion, as 

PS,a,a+i(e,e',e") = PS(e)Pi^)(e')^i+\(e") 

+SP [PS(e)gHi(e', e") + Pii_\(e")g?ii,„(e, e'); , 



Plugging this expression into equation (]34l) . we obtain an equation containing the terms 
of equation (l36l) . but this time with additional terms involving g's. Though this equation 
is complicated and yields a priori no simple explicit solution for q, it is consistent with 
the fact that j does not have the gradient property. 

Nonetheless, finding the exact expression of q is not necessary for our sake as 
the symmetry of q, equation (H3l) . turns out to be enough information on the two-cell 
distribution function to compute the stationary heat current fl27|) . Indeed, we may write 
the stationary current as the average 

Jh = / dxd|/e-(^+^) |l - ^[y - q{x,y)]\j . (45) 
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Notice from equation ([H]) that j{x/P,y/P) = P~^^'^j{x,y). Thus, to 0{SP), there are 
only two possible contributions to the heat current that need to be taken into account. 
Namely, 

dxdye~^^~^'^^yj{x, y) 

(46) 

dxdye ^^~^^^q{x,y)j{x,y) . 



However the second of these expressions identically vanishes since j is an ant i- symmetric 
function of its arguments and q, as we saw in equation ( H3l) . is symmetric. The steady 
state current thus becomes 

"^H = I dxdye-^^+y^yj{x, y) , (47) 



^5/2 

which, by symmetry, can be rewritten as 

"^H = I dxd?/e~(^+^)(x - y)jix,y) 



dxdye ^^~^y\x - y^i^ix, y) 



2/55/2 
1 5(3 



^5/2 ' 

= -VtST. (48) 

where we consecutively used equations (El) in the second line, (fTOi) in the third one, and 
the identity 6(3 = —ST/T"^ in the last one. Equation (HHj) is nothing but Fourier's law 
(j26ll . with conductivity 

K = Vf, (49) 

confirming equation fl23l) . Thus, even though there are 0{S(3) contributions to the 
stationary state which arise from the two-point distribution function, the only actual 
contribution to heat transport arises from the average with respect to the local 
equilibrium part. This result also confirms our claim that the cross-correlations on the 
RHS of equation vanish in the large system limit and that only the auto-correlation 
part (|22|) contributes to the heat conductivity. 

Before turning to numerical results in the next section, we should note that our 
arguments though compelling do not constitute a formal proof that the distributions (!38ll 
with q symmetric are the unique stationary solutions of equations ( l25l) . ( 1341) consistent 
with the non-equilibrium boundary conditions. This is in fact a much deeper problem 
that goes beyond the scope of the present paper. What we have shown rather is that, 
given that the system is in contact with thermal baths at distinct temperatures, there 
exist stationary states of the form (138|1 with associated heat current obeying Fourier's 
law (HHl) with the corresponding conductivity (H9l) . It is a natural conjecture to assume 
that, given the baths temperatures, the stationary state fl38l) is unique. We turn to 
numerical simulations to assess that claim. 
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7. Numerical computation 



Equations (!23l) and ( 1491) can be verified by direct numerical computations of tlie master 
equation ([1]). Tliis is acliieved by adapting Gillespie's algorithm [6], [7], originally designed 
to simulate chemical or biochemical systems of reactions, to our stochastic equation. 
The Monte-Carlo step here necessitates two random trials. The first random number 
determines the time that will elapse until the next energy exchange event. The second 
one determines which one out of all the possible pairs of cells will perform an exchange 
of energy and how much energy will be exchanged between them. 

Thus, given a configuration {ei, . . . , cn} of energies at each one of the N cells of the 
system, and, whether in the presence of periodic boundary conditions (PBC), used to 
simulate an isolated system, or thermal boundary conditions (TBC), used to simulate 
a system with both ends in contact with thermostats at distinct temperatures T_ and 
T+, and, in that case, given a pre-specified bath relaxation rate i^Bath associated with 
the thermal baths, the time to the next collision is a random number with Poisson 
distribution whose relaxation rate is specified by the sum of all collision frequencies, 

7V-1 

7=^2 ^(^ri, e„+i) + z/(ejv, d) , (PBC) , 

(50) 

7 = i^Bath ( \/t2 + + ^('^ri, e„+i) , (TBC) . 

^ ^ n=0 

Given the rate 7, one draws a second random number x uniformly distributed on (0, 7), 
and determines the pair (n, n + 1) of cells which interact according to 

n— 1 71 

J2 '^(ei, ei+i) <X<Y1 ^(^i^ ^i+i) ' (PBC) , 

^=l ^=l ^5^^ 
n— 1 71 ^ ' 

'^Bath + '^i^i^ <X < '^BathyT^ + ^ u{ei, e^+i) , (TBC) . 

1=0 i=0 

In the latter equation, one understands that x < ^BathVTZ means updating the energy 
of the left bath, also associated with i = 0, and, similarly, x > ^Bathy/TZ+YliLo ^i^i^ ^i+i) 
means updating the energy of the right bath, associated with cell + 1. How much 
energy is exchanged between the cells n and n + 1 is then determined by finding the 
energy 77 which solves the equation 

71—1 

-en 

Here we assumed PBC. The transposition to TBC is immediate. For the sake of solving 
this equation, going back to the expression of the collision frequency (Il])-([5]), we compute 
the partial integrals 



1 n 

/ dri'W{ea, eb\ea - r/', + r]') 

jli ^b) "'-et 
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'-K\ 



2El ^]-ibzl2.K{ ^ 
I \'b 



<^b 2E[ is. 

^6 



-K[ ^ 

^b 



2E[ ^ 



"b V 



-eb < r] < ea - eb . 
ea - eb < T] < , 
0<r]<ea, 



(53) 



and use a root finder routine -in our case the subroutine rtf Isp in [9]- in order to find 
the solution rj to equation fl52l) . 

We notice that the elhptic K{x) functions which appear in equation (l53l) . are 
multiphed by 1 — x, which takes care of the logarithmic divergence of at a; = 1, 



hm 



Kix) 



1 , 16 
- log 

2 ^ 1 -a; 



0. 



(54) 



7.1. Periodic Boundary Conditions 



As a simple test of the algorithm, we compute the equilibrium average of the collision 
frequency and compare it to the results of a numerical integration of equation ([9]) for 
different system sizes and energy per particle taken to be unity. The results, displayed 
in the top panel of figure [2|, are in excellent agreement. 

Consider then the Helfand moment ( fT6l) . We present in the bottom panel of figure 
[2] the results of numerical computations of equation ( |T9i) . obtained for different sizes 
N of the system, using the numerical scheme described above. As can be seen from 
the figure, lim^^oo l/(2A^)(Aif(r„)^/r„), has C(l/A^) corrections to z/b which are due 
to cross-correlations of the form (//(e^., efc.+i)?7(efc^. , ejt^+i)) that persist for finite N. The 
infinite N extrapolation of our data is obtained by linear regression and yields 

k/z/b = 0.997 ± 0.004 , (55) 

in precise agreement with equation (123!) and our argument that only static correlations 
contribute to the thermal conductivity. 



7.2. Thermal Boundary Conditions 

We simulate thermal boundary conditions according to the scheme described above, 
with bath relaxation frequency i^Bath = 100, fixing the baths temperatures to T„ = 0.5 
and T+ = 1.5. The heat conductivity is then evaluated according to equation (!26ll 
by computing the average heat flux, divided by the temperature gradient, which is 
decreased by increasing the system size N . The results, displayed on the left panel of 
figure [HI yield the ratio 

k/z/b = 1.0002 ±3 10^'S (56) 

which provides a very nice confirmation of equation (jUj). 
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Figure 2. (Top) Collision frequency computed under equilibrium conditions with 
energy E ^ N and for different system sizes N, ranging from = 3 to = 100, 
compared to the results of a numerical integration of ([9|) (solid red line). The N ^ oo 
asymptotic value is 1 according to (H]). The agreement is spectacular. (Bottom) RHS 
of equation (fTQ]) . computed for the same parameters set as above. The ratio k/i^b is 
the infinite N limit of these data, which can be evaluated by linear regression (solid 
red line). Fitting all the data points corresponding to iV > 10 and weighting the data 
according to the sizes of their error bars, we obtain k/i^b = 0.997 ± 0.004. 
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Figure 3. (Left) Ratio between the heat current ((27|) and local temperature gradient 
under thermal boundary conditions at temperatures r_|_ — 1.5 and T_ = 0.5, computed 
for different system sizes N, ranging from = 1 to = 30. The ratio k/vb 
is computed from the infinite N limit of these data, evaluated by linear regression 
(solid red line). Fitting the data points weighed according to the sizes of their error 
bars, we obtain k/vb — 1.0002 ± 310^'*. (Right) Corresponding temperature profiles 
with stochastically thermalised cells at n = (r_) and n = N + 1 (r+). The thick 
black line corresponds to profile (l57|) expected from Fourier's law. The data is barely 
distinguishable from this curve when N is sufficiently large. 



The corresponding temperature profiles which, according to Fourier's law, are 
expected to be 



T 



n 



{t: 



,3/2 



3/2n 



1 2/3 



(57) 



are shown on the right panel of figure [3l 



8. Conclusions 



To summarize, we have obtained a systematic derivation of Fourier's law and computed 
the heat conductivity of a class of stochastic systems describing, at a mesoscopic level, 
the rare energy exchanges in Hamiltonian systems of two-dimensional confined particles 
in interaction that were introduced in [Tl|2]. 

A remarkable feature of our approach is that the structure of the stationary state 
is determined by the geometry of the system and not by the actual form of the kernel 
which specifies the nature of the interactions between the system cells. Indeed our 
derivation of the two-cell distribution function only used the nearest neighbour nature 
of the interaction. This property alone justifies that the leading part of the distribution 
function, meaning including linear order in the local temperature gradients, is the local 
thermal equilibrium plus a two-cell function, symmetric with respect to its arguments. 
What precise form this function has depends on the actual process under consideration. 
Whatever this form, its symmetry suffices to justify that the only contributions to the 
heat current come from the local thermal equilibrium part, at least to the extent that 
the current is computed to linear order in the local temperature gradients. 
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We therefore conclude the same property, namely that the heat conductivity can be 
obtained through averages with respect to the local thermal equilibrium distributions, 
should hold for systems similarly described by an equation of the form ([T]), sustaining 
non-equilibrium stationary states with symmetric corrections of the form (!38l) . 

In particular, we observe that, whether or not such systems obey the gradient 
condition, the transport coefficient is given from the diffusion of Helfand moment, or 
equivalently the Green-Kubo formula, through a static average only. 

Let us mention that, as observed in [10], the class of mechanical systems with 
confined particles in interaction that share the transport properties of the systems 
considered in [U [2] is larger than the class of semi- dispersing billiards envisaged there. 
In particular, the systems studied in [10] do not a priori obey a master equation of the 
form ([T]). We expect that the BBGKY hierarchy applied to the phase-space distributions 
of these systems will have solutions with properties similar to the solutions found here, 
namely a local thermal equilibrium part plus a symmetric correction which does not 
contribute to the heat current. 
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